Near Analysis Models



What Was Tested: Do crop yields increase from GCOs globally as a function of fertilizer? Note that since some crops have the same GCOs their distance bands are calculated once.



What This Allows Me To Do: Identify if yield-distance relationships (near analysis) are a function of evolutionary escape or other fertilizer.



What Are The Model Assumptions: GLS models are used as they have less assumptions about error structure and heteroscadicity.

library(rgdal)
library(maptools)
library(rgeos)
library(ggplot2)
library(dplyr)
library(reshape2)

#ggplot graphic parameters
theme_justin<-theme_bw() +theme(axis.line = element_line(colour = "black"),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.background = element_blank())

#Import Fishnet
fishnet<- readOGR(dsn= "/Users/justinstewart/Dropbox/Collaborations/TobyKiers/Crop_Productivity_GCO/Analysis/GIS/Near_Fertilizer/", layer="Near_Fertilizer")
## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/justinstewart/Dropbox/Collaborations/TobyKiers/Crop_Productivity_GCO/Analysis/GIS/Near_Fertilizer", layer: "Near_Fertilizer"
## with 37445 features
## It has 16 fields
## Integer64 fields read as strings:  OBJECTID
### Get Lat-Long Baby Girl and apply metadata
xy <- coordinates(fishnet)
metadata <- fishnet@data
sapply(metadata, "class")
##    OBJECTID     COUNTRY  Fishnet_ID  Fishnet__1  barley_Fer  cassava_Fe 
## "character" "character"   "integer"   "numeric"   "numeric"   "numeric" 
##  groundnut_  maize_Fert  millet_Fer  rapeseed_F  rice_FertM  rye_FertMe 
##   "numeric"   "numeric"   "numeric"   "numeric"   "numeric"   "numeric" 
##  sorghum_Fe  soybean_Fe  sunflower_  wheat_Fert 
##   "numeric"   "numeric"   "numeric"   "numeric"
plot(fishnet) #Fishnet only contains pixels with fertilizer input

## Barley

barley<- readOGR(dsn= "/Users/justinstewart/Dropbox/Collaborations/TobyKiers/Crop_Productivity_GCO/Analysis/GIS/Near_Fertilizer/", layer="Barley_Wheat_Centroid") #Import GCO Centroid
## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/justinstewart/Dropbox/Collaborations/TobyKiers/Crop_Productivity_GCO/Analysis/GIS/Near_Fertilizer", layer: "Barley_Wheat_Centroid"
## with 1 features
## It has 2 fields
## Integer64 fields read as strings:  OBJECTID ORIG_FID
barley.c<-coordinates(barley) #Extract Lat-Long
distances <- spDistsN1(xy,barley.c, longlat = FALSE) #Calculate Near Distances
data<-cbind(data.frame(metadata),data.frame(distances)) #Merge into a DF
data$rescale_ND<-data$distances/(1000*1000) #Rescale distances 
data<-data %>% filter(barley_Fer > 0, na.rm=TRUE) #select Fertilizer > 0 

#Pivot Table To ID Top Fertilizer Applications by Country Sum 
barley.pivot<-dcast(data, COUNTRY ~., value.var="barley_Fer", fun.aggregate=sum)
barley.pivot<-arrange(barley.pivot,.)
tail(barley.pivot, 3) #ID Top 3 Fertilizer Countries
##           COUNTRY        .
## 161        Canada 20686.84
## 162 United States 22053.82
## 163         China 52384.10
sum(barley.pivot$.) #Total Fertilizer
## [1] 299702.4
#Filter Countries
china<-filter(data, COUNTRY == "China") 
china.Perc<-sum(china$barley_Fer)/sum(barley.pivot$.) 
print(china.Perc) # % of Global Fertilizer For Crop 
## [1] 0.174787
china.color<-"#1665AF"
summary(china)
##    OBJECTID           COUNTRY            Fishnet_ID      Fishnet__1   
##  Length:1417        Length:1417        Min.   :  464   Min.   :  464  
##  Class :character   Class :character   1st Qu.:48702   1st Qu.:48702  
##  Mode  :character   Mode  :character   Median :57638   Median :57638  
##                                        Mean   :54982   Mean   :54982  
##                                        3rd Qu.:64219   3rd Qu.:64219  
##                                        Max.   :76018   Max.   :76018  
##    barley_Fer        cassava_Fe        groundnut_         maize_Fert      
##  Min.   : 0.3023   Min.   : 0.5672   Min.   :  0.5672   Min.   :  0.5672  
##  1st Qu.:33.6040   1st Qu.:36.2767   1st Qu.: 30.5924   1st Qu.: 52.4740  
##  Median :35.4694   Median :37.5641   Median : 40.0470   Median : 66.0109  
##  Mean   :36.9683   Mean   :39.1958   Mean   : 40.7972   Mean   : 61.0785  
##  3rd Qu.:44.4152   3rd Qu.:47.9056   3rd Qu.: 52.7165   3rd Qu.: 75.1105  
##  Max.   :92.0236   Max.   :62.6335   Max.   :126.9427   Max.   :227.1572  
##    millet_Fer         rapeseed_F         rice_FertM         rye_FertMe     
##  Min.   :  0.5672   Min.   :  0.5672   Min.   :  0.4882   Min.   : 0.5672  
##  1st Qu.:  8.0391   1st Qu.: 54.1873   1st Qu.: 26.4062   1st Qu.:31.9068  
##  Median : 33.5524   Median : 58.4773   Median : 62.7289   Median :35.4209  
##  Mean   : 30.2364   Mean   : 60.2659   Mean   : 56.8110   Mean   :36.6549  
##  3rd Qu.: 44.1182   3rd Qu.: 71.7545   3rd Qu.: 79.1286   3rd Qu.:42.6427  
##  Max.   :105.8587   Max.   :151.1022   Max.   :158.1863   Max.   :92.6647  
##    sorghum_Fe         soybean_Fe         sunflower_        wheat_Fert      
##  Min.   :  0.5672   Min.   :  0.5672   Min.   : 0.5672   Min.   :  0.5481  
##  1st Qu.: 39.7837   1st Qu.: 30.7152   1st Qu.:18.3559   1st Qu.: 52.8928  
##  Median : 51.7188   Median : 40.5146   Median :37.4452   Median : 59.8756  
##  Mean   : 48.0861   Mean   : 40.9564   Mean   :33.6039   Mean   : 61.0944  
##  3rd Qu.: 66.6392   3rd Qu.: 49.7050   3rd Qu.:47.3547   3rd Qu.: 71.2560  
##  Max.   :116.3278   Max.   :111.8416   Max.   :89.3264   Max.   :148.6335  
##    distances          rescale_ND    
##  Min.   : 9387729   Min.   : 9.388  
##  1st Qu.:11428619   1st Qu.:11.429  
##  Median :12695259   Median :12.695  
##  Mean   :12671117   Mean   :12.671  
##  3rd Qu.:13863752   3rd Qu.:13.864  
##  Max.   :15943817   Max.   :15.944
usa<-filter(data, COUNTRY == "United States")
usa.Perc<-sum(usa$barley_Fer)/sum(barley.pivot$.) 
print(usa.Perc) # % of Global Fertilizer For Crop 
## [1] 0.07358571
usa.color<-"#8F4A33"

canada<-filter(data, COUNTRY == "Canada")
canada.Perc<-sum(canada$barley_Fer)/sum(barley.pivot$.) 
print(canada.Perc) # % of Global Fertilizer For Crop 
## [1] 0.06902458
canada.color<-"#98C24B"

#Plot 
ggplot(data, aes(x=rescale_ND,y=log10(data$barley_Fer))) + 
      geom_point(alpha=0.3) +theme_justin + 
      geom_point(data=china, aes(x=rescale_ND, y=log10(china$barley_Fer)), color=china.color, alpha=0.1) + 
      geom_point(data=usa, aes(x=rescale_ND, y=log10(usa$barley_Fer)), color=usa.color, alpha=0.2) + 
      geom_point(data=canada, aes(x=rescale_ND, y=log10(canada$barley_Fer)), color=canada.color, alpha=0.1) + 
      geom_smooth(method="lm", color="red") +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Fertilizer Input")

###  Model - Barley Fertilizer Near
library(nlme)
data$logBarley<-log10(data$barley_Fer)
# regular OLS no variance structure
barley.ols <- gls(logBarley ~ rescale_ND, data = data)
# varFixed (variance changes linearly with X)
barley.fixed <- update(barley.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
barley.power <- update(barley.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
barley.exp <- update(barley.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
barley.ConstPower <- update(barley.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(barley.ols, barley.fixed, barley.power, barley.ConstPower, barley.exp)
##                   df      AIC
## barley.ols         3 29265.90
## barley.fixed       3 30503.28
## barley.power       4 29238.88
## barley.ConstPower  5 29240.88
## barley.exp         4 29261.54
summary(barley.power)
## Generalized least squares fit by REML
##   Model: logBarley ~ rescale_ND 
##   Data: data 
##        AIC      BIC    logLik
##   29238.88 29269.56 -14615.44
## 
## Variance function:
##  Structure: Power of variance covariate
##  Formula: ~rescale_ND 
##  Parameter estimates:
##      power 
## 0.06477432 
## 
## Coefficients:
##                 Value   Std.Error  t-value p-value
## (Intercept) 0.6228270 0.011735326 53.07284       0
## rescale_ND  0.0321897 0.001223496 26.30962       0
## 
##  Correlation: 
##            (Intr)
## rescale_ND -0.911
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -2.9688596 -0.8093478  0.1257739  0.7259930  3.7542355 
## 
## Residual standard error: 0.5323013 
## Degrees of freedom: 15830 total; 15828 residual

Cassava

cassava<- readOGR(dsn= "/Users/justinstewart/Dropbox/Collaborations/TobyKiers/Crop_Productivity_GCO/Analysis/GIS/Near_Fertilizer/", layer="Cassava_Groundnut_Centroid") #Import GCO Centroid
## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/justinstewart/Dropbox/Collaborations/TobyKiers/Crop_Productivity_GCO/Analysis/GIS/Near_Fertilizer", layer: "Cassava_Groundnut_Centroid"
## with 1 features
## It has 2 fields
## Integer64 fields read as strings:  OBJECTID ORIG_FID
cassava.c<-coordinates(cassava) #Extract Lat-Long
distances <- spDistsN1(xy,cassava.c, longlat = FALSE) #Calculate Near Distances
data<-cbind(data.frame(metadata),data.frame(distances)) #Merge into a DF
data$rescale_ND<-data$distances/(1000*1000) #Rescale distances 
data<-data %>% filter(cassava_Fe > 0, na.rm=TRUE) #select Fertilizer > 0 

#Pivot Table To ID Top Fertilizer Applications by Country Sum 
cassava.pivot<-dcast(data, COUNTRY ~., value.var="cassava_Fe", fun.aggregate=sum)
cassava.pivot<-arrange(cassava.pivot,.)
tail(cassava.pivot, 3) #ID Top 3 Fertilizer Countries
##           COUNTRY         .
## 161 United States  5867.965
## 162        Brazil 28021.176
## 163         China 55540.387
sum(cassava.pivot$.) #Total Fertilizer
## [1] 154832.8
#Filter Countries
china<-filter(data, COUNTRY == "China") 
china.Perc<-sum(china$cassava_Fer)/sum(cassava.pivot$.) 
print(china.Perc) # % of Global Fertilizer For Crop 
## [1] 0
china.color<-"#1665AF"
summary(china)
##    OBJECTID           COUNTRY            Fishnet_ID      Fishnet__1   
##  Length:1417        Length:1417        Min.   :  464   Min.   :  464  
##  Class :character   Class :character   1st Qu.:48702   1st Qu.:48702  
##  Mode  :character   Mode  :character   Median :57638   Median :57638  
##                                        Mean   :54982   Mean   :54982  
##                                        3rd Qu.:64219   3rd Qu.:64219  
##                                        Max.   :76018   Max.   :76018  
##    barley_Fer        cassava_Fe        groundnut_         maize_Fert      
##  Min.   : 0.3023   Min.   : 0.5672   Min.   :  0.5672   Min.   :  0.5672  
##  1st Qu.:33.6040   1st Qu.:36.2767   1st Qu.: 30.5924   1st Qu.: 52.4740  
##  Median :35.4694   Median :37.5641   Median : 40.0470   Median : 66.0109  
##  Mean   :36.9683   Mean   :39.1958   Mean   : 40.7972   Mean   : 61.0785  
##  3rd Qu.:44.4152   3rd Qu.:47.9056   3rd Qu.: 52.7165   3rd Qu.: 75.1105  
##  Max.   :92.0236   Max.   :62.6335   Max.   :126.9427   Max.   :227.1572  
##    millet_Fer         rapeseed_F         rice_FertM         rye_FertMe     
##  Min.   :  0.5672   Min.   :  0.5672   Min.   :  0.4882   Min.   : 0.5672  
##  1st Qu.:  8.0391   1st Qu.: 54.1873   1st Qu.: 26.4062   1st Qu.:31.9068  
##  Median : 33.5524   Median : 58.4773   Median : 62.7289   Median :35.4209  
##  Mean   : 30.2364   Mean   : 60.2659   Mean   : 56.8110   Mean   :36.6549  
##  3rd Qu.: 44.1182   3rd Qu.: 71.7545   3rd Qu.: 79.1286   3rd Qu.:42.6427  
##  Max.   :105.8587   Max.   :151.1022   Max.   :158.1863   Max.   :92.6647  
##    sorghum_Fe         soybean_Fe         sunflower_        wheat_Fert      
##  Min.   :  0.5672   Min.   :  0.5672   Min.   : 0.5672   Min.   :  0.5481  
##  1st Qu.: 39.7837   1st Qu.: 30.7152   1st Qu.:18.3559   1st Qu.: 52.8928  
##  Median : 51.7188   Median : 40.5146   Median :37.4452   Median : 59.8756  
##  Mean   : 48.0861   Mean   : 40.9564   Mean   :33.6039   Mean   : 61.0944  
##  3rd Qu.: 66.6392   3rd Qu.: 49.7050   3rd Qu.:47.3547   3rd Qu.: 71.2560  
##  Max.   :116.3278   Max.   :111.8416   Max.   :89.3264   Max.   :148.6335  
##    distances          rescale_ND    
##  Min.   : 9387843   Min.   : 9.388  
##  1st Qu.:11428732   1st Qu.:11.429  
##  Median :12695373   Median :12.695  
##  Mean   :12671230   Mean   :12.671  
##  3rd Qu.:13863866   3rd Qu.:13.864  
##  Max.   :15943931   Max.   :15.944
brazil<-filter(data, COUNTRY == "Brazil")
brazil.Perc<-sum(brazil$cassava_Fer)/sum(cassava.pivot$.) 
print(brazil.Perc) # % of Global Fertilizer For Crop 
## [1] 0
brazil.color<-"#E28D15"

usa<-filter(data, COUNTRY == "United States")
usa.Perc<-sum(usa$cassava_Fer)/sum(cassava.pivot$.) 
print(usa.Perc) # % of Global Fertilizer For Crop 
## [1] 0
usa.color<-"#8F4A33"


#Plot 
ggplot(data, aes(x=rescale_ND,y=log10(data$cassava_Fe))) + 
      geom_point(alpha=0.3) +theme_justin + 
      geom_point(data=china, aes(x=rescale_ND, y=log10(china$cassava_Fe)), color=china.color, alpha=0.1) + 
      geom_point(data=brazil, aes(x=rescale_ND, y=log10(brazil$cassava_Fe)), color=brazil.color, alpha=0.2) + 
      geom_point(data=usa, aes(x=rescale_ND, y=log10(usa$cassava_Fe)), color=usa.color, alpha=0.1) + 
      geom_smooth(method="lm", color="red") +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Fertilizer Input")

###  Model - Cassava Fertilizer Near
library(nlme)
data$logCassava<-log10(data$cassava_Fe)
# regular OLS no variance structure
cassava.ols <- gls(logCassava ~ rescale_ND, data = data)
# varFixed (variance changes linearly with X)
cassava.fixed <- update(cassava.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
cassava.power <- update(cassava.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
cassava.exp <- update(cassava.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
cassava.ConstPower <- update(cassava.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(cassava.ols, cassava.fixed, cassava.power, cassava.ConstPower, cassava.exp)
##                    df      AIC
## cassava.ols         3 27154.72
## cassava.fixed       3 28539.68
## cassava.power       4 26874.13
## cassava.ConstPower  5 26638.64
## cassava.exp         4 26705.41
summary(cassava.exp)
## Generalized least squares fit by REML
##   Model: logCassava ~ rescale_ND 
##   Data: data 
##        AIC      BIC    logLik
##   26705.41 26736.09 -13348.71
## 
## Variance function:
##  Structure: Exponential of variance covariate
##  Formula: ~rescale_ND 
##  Parameter estimates:
##      expon 
## 0.02975004 
## 
## Coefficients:
##                  Value   Std.Error  t-value p-value
## (Intercept)  0.6720116 0.010301379 65.23511       0
## rescale_ND  -0.0070547 0.001147294 -6.14896       0
## 
##  Correlation: 
##            (Intr)
## rescale_ND -0.904
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -2.7679892 -0.7520635 -0.2185703  0.6662134  4.0932392 
## 
## Residual standard error: 0.4298106 
## Degrees of freedom: 15830 total; 15828 residual

Groundnut

groundnut<- readOGR(dsn= "/Users/justinstewart/Dropbox/Collaborations/TobyKiers/Crop_Productivity_GCO/Analysis/GIS/Near_Fertilizer/", layer="Cassava_Groundnut_Centroid") #Import GCO Centroid
## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/justinstewart/Dropbox/Collaborations/TobyKiers/Crop_Productivity_GCO/Analysis/GIS/Near_Fertilizer", layer: "Cassava_Groundnut_Centroid"
## with 1 features
## It has 2 fields
## Integer64 fields read as strings:  OBJECTID ORIG_FID
groundnut.c<-coordinates(groundnut) #Extract Lat-Long
distances <- spDistsN1(xy,groundnut.c, longlat = FALSE) #Calculate Near Distances
data<-cbind(data.frame(metadata),data.frame(distances)) #Merge into a DF
data$rescale_ND<-data$distances/(1000*1000) #Rescale distances 
data<-data %>% filter(groundnut_> 0, na.rm=TRUE) #select Fertilizer > 0 

#Pivot Table To ID Top Fertilizer Applications by Country Sum 
groundnut.pivot<-dcast(data, COUNTRY ~., value.var="groundnut_", fun.aggregate=sum)
groundnut.pivot<-arrange(groundnut.pivot,.)
tail(groundnut.pivot, 3) #ID Top 3 Fertilizer Countries
##           COUNTRY        .
## 161 United States 10164.94
## 162         India 10557.89
## 163         China 57809.56
sum(groundnut.pivot$.) #Total Fertilizer
## [1] 162718.7
#Filter Countries
china<-filter(data, COUNTRY == "China") 
china.Perc<-sum(china$groundnut_Fer)/sum(groundnut.pivot$.) 
print(china.Perc) # % of Global Fertilizer For Crop 
## [1] 0
china.color<-"#1665AF"
summary(china)
##    OBJECTID           COUNTRY            Fishnet_ID      Fishnet__1   
##  Length:1417        Length:1417        Min.   :  464   Min.   :  464  
##  Class :character   Class :character   1st Qu.:48702   1st Qu.:48702  
##  Mode  :character   Mode  :character   Median :57638   Median :57638  
##                                        Mean   :54982   Mean   :54982  
##                                        3rd Qu.:64219   3rd Qu.:64219  
##                                        Max.   :76018   Max.   :76018  
##    barley_Fer        cassava_Fe        groundnut_         maize_Fert      
##  Min.   : 0.3023   Min.   : 0.5672   Min.   :  0.5672   Min.   :  0.5672  
##  1st Qu.:33.6040   1st Qu.:36.2767   1st Qu.: 30.5924   1st Qu.: 52.4740  
##  Median :35.4694   Median :37.5641   Median : 40.0470   Median : 66.0109  
##  Mean   :36.9683   Mean   :39.1958   Mean   : 40.7972   Mean   : 61.0785  
##  3rd Qu.:44.4152   3rd Qu.:47.9056   3rd Qu.: 52.7165   3rd Qu.: 75.1105  
##  Max.   :92.0236   Max.   :62.6335   Max.   :126.9427   Max.   :227.1572  
##    millet_Fer         rapeseed_F         rice_FertM         rye_FertMe     
##  Min.   :  0.5672   Min.   :  0.5672   Min.   :  0.4882   Min.   : 0.5672  
##  1st Qu.:  8.0391   1st Qu.: 54.1873   1st Qu.: 26.4062   1st Qu.:31.9068  
##  Median : 33.5524   Median : 58.4773   Median : 62.7289   Median :35.4209  
##  Mean   : 30.2364   Mean   : 60.2659   Mean   : 56.8110   Mean   :36.6549  
##  3rd Qu.: 44.1182   3rd Qu.: 71.7545   3rd Qu.: 79.1286   3rd Qu.:42.6427  
##  Max.   :105.8587   Max.   :151.1022   Max.   :158.1863   Max.   :92.6647  
##    sorghum_Fe         soybean_Fe         sunflower_        wheat_Fert      
##  Min.   :  0.5672   Min.   :  0.5672   Min.   : 0.5672   Min.   :  0.5481  
##  1st Qu.: 39.7837   1st Qu.: 30.7152   1st Qu.:18.3559   1st Qu.: 52.8928  
##  Median : 51.7188   Median : 40.5146   Median :37.4452   Median : 59.8756  
##  Mean   : 48.0861   Mean   : 40.9564   Mean   :33.6039   Mean   : 61.0944  
##  3rd Qu.: 66.6392   3rd Qu.: 49.7050   3rd Qu.:47.3547   3rd Qu.: 71.2560  
##  Max.   :116.3278   Max.   :111.8416   Max.   :89.3264   Max.   :148.6335  
##    distances          rescale_ND    
##  Min.   : 9387843   Min.   : 9.388  
##  1st Qu.:11428732   1st Qu.:11.429  
##  Median :12695373   Median :12.695  
##  Mean   :12671230   Mean   :12.671  
##  3rd Qu.:13863866   3rd Qu.:13.864  
##  Max.   :15943931   Max.   :15.944
india<-filter(data, COUNTRY == "India")
india.Perc<-sum(canada$groundnut_Fer)/sum(groundnut.pivot$.) 
print(india.Perc) # % of Global Fertilizer For Crop 
## [1] 0
india.color<-"#078D40"

usa<-filter(data, COUNTRY == "United States")
usa.Perc<-sum(usa$groundnut_Fer)/sum(groundnut.pivot$.) 
print(usa.Perc) # % of Global Fertilizer For Crop 
## [1] 0
usa.color<-"#8F4A33"

#Plot 
ggplot(data, aes(x=rescale_ND,y=log10(data$groundnut_))) + 
      geom_point(alpha=0.3) +theme_justin + 
      geom_point(data=china, aes(x=rescale_ND, y=log10(china$groundnut_)), color=china.color, alpha=0.1) + 
      geom_point(data=india, aes(x=rescale_ND, y=log10(india$groundnut_)), color=usa.color, alpha=0.2) + 
      geom_point(data=usa, aes(x=rescale_ND, y=log10(usa$groundnut_)), color=canada.color, alpha=0.1) + 
      geom_smooth(method="lm", color="red") +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Fertilizer Input")

###  Model - Groundnut Fertilizer Near
library(nlme)
data$logGroundnut<-log10(data$groundnut_)
# regular OLS no variance structure
groundnut.ols <- gls(logGroundnut ~ rescale_ND, data = data)
# varFixed (variance changes linearly with X)
groundnut.fixed <- update(groundnut.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
groundnut.power <- update(groundnut.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
groundnut.exp <- update(groundnut.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
groundnut.ConstPower <- update(groundnut.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(groundnut.ols, groundnut.fixed, groundnut.power, groundnut.ConstPower, groundnut.exp)
##                      df      AIC
## groundnut.ols         3 27883.16
## groundnut.fixed       3 28029.44
## groundnut.power       4 27182.17
## groundnut.ConstPower  5 26837.89
## groundnut.exp         4 26920.76
summary(groundnut.exp)
## Generalized least squares fit by REML
##   Model: logGroundnut ~ rescale_ND 
##   Data: data 
##        AIC      BIC    logLik
##   26920.76 26951.44 -13456.38
## 
## Variance function:
##  Structure: Exponential of variance covariate
##  Formula: ~rescale_ND 
##  Parameter estimates:
##      expon 
## 0.04343479 
## 
## Coefficients:
##                  Value   Std.Error  t-value p-value
## (Intercept)  0.6592527 0.009925003 66.42343  0.0000
## rescale_ND  -0.0031826 0.001154617 -2.75644  0.0059
## 
##  Correlation: 
##            (Intr)
## rescale_ND -0.898
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -2.99346011 -0.78095060 -0.06551759  0.73297285  4.28379085 
## 
## Residual standard error: 0.38253 
## Degrees of freedom: 15830 total; 15828 residual

Sorghum

sorghum<- readOGR(dsn= "/Users/justinstewart/Dropbox/Collaborations/TobyKiers/Crop_Productivity_GCO/Analysis/GIS/Near_Fertilizer/", layer="Sorghum_Centroid") #Import GCO Centroid
## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/justinstewart/Dropbox/Collaborations/TobyKiers/Crop_Productivity_GCO/Analysis/GIS/Near_Fertilizer", layer: "Sorghum_Centroid"
## with 1 features
## It has 2 fields
## Integer64 fields read as strings:  OBJECTID ORIG_FID
sorghum.c<-coordinates(sorghum) #Extract Lat-Long
distances <- spDistsN1(xy,sorghum.c, longlat = FALSE) #Calculate Near Distances
data<-cbind(data.frame(metadata),data.frame(distances)) #Merge into a DF
data$rescale_ND<-data$distances/(1000*1000) #Rescale distances 
data<-data %>% filter(sorghum_Fe > 0, na.rm=TRUE) #select Fertilizer > 0 
summary(data)
##    OBJECTID           COUNTRY            Fishnet_ID      Fishnet__1   
##  Length:15830       Length:15830       Min.   :   19   Min.   :   19  
##  Class :character   Class :character   1st Qu.:45702   1st Qu.:45702  
##  Mode  :character   Mode  :character   Median :53450   Median :53450  
##                                        Mean   :52349   Mean   :52349  
##                                        3rd Qu.:62261   3rd Qu.:62261  
##                                        Max.   :80977   Max.   :80977  
##    barley_Fer          cassava_Fe         groundnut_         maize_Fert      
##  Min.   :   0.1264   Min.   :  0.1027   Min.   :  0.1027   Min.   :  0.1264  
##  1st Qu.:   2.3293   1st Qu.:  1.4786   1st Qu.:  1.4558   1st Qu.:  4.8776  
##  Median :  10.3748   Median :  3.1312   Median :  3.9874   Median : 15.3926  
##  Mean   :  18.9326   Mean   :  9.7810   Mean   : 10.2791   Mean   : 29.1414  
##  3rd Qu.:  27.3848   3rd Qu.:  9.6465   3rd Qu.: 10.8990   3rd Qu.: 45.3560  
##  Max.   :1255.3557   Max.   :323.1332   Max.   :261.9184   Max.   :821.1623  
##    millet_Fer         rapeseed_F         rice_FertM         rye_FertMe      
##  Min.   :  0.1027   Min.   :  0.1264   Min.   :  0.1027   Min.   :  0.1877  
##  1st Qu.:  1.4394   1st Qu.:  2.1127   1st Qu.:  1.5730   1st Qu.:  1.9372  
##  Median :  3.0850   Median :  5.8253   Median :  4.5741   Median :  4.4323  
##  Mean   :  8.8622   Mean   : 20.9648   Mean   : 16.4947   Mean   : 12.5756  
##  3rd Qu.:  9.4221   3rd Qu.: 31.7499   3rd Qu.: 22.4244   3rd Qu.: 16.7698  
##  Max.   :278.4499   Max.   :166.8706   Max.   :487.8604   Max.   :158.1257  
##    sorghum_Fe         soybean_Fe         sunflower_         wheat_Fert      
##  Min.   :  0.1027   Min.   :  0.1027   Min.   :  0.1087   Min.   :  0.1264  
##  1st Qu.:  1.6315   1st Qu.:  1.8021   1st Qu.:  1.7595   1st Qu.:  2.6795  
##  Median :  5.5436   Median :  3.9118   Median :  4.5399   Median : 14.3677  
##  Mean   : 13.7190   Mean   : 12.1588   Mean   : 12.0596   Mean   : 24.6181  
##  3rd Qu.: 17.7258   3rd Qu.: 17.2672   3rd Qu.: 13.1712   3rd Qu.: 39.8120  
##  Max.   :955.8894   Max.   :139.3965   Max.   :542.2973   Max.   :530.9011  
##    distances          rescale_ND     
##  Min.   :  597524   Min.   : 0.5975  
##  1st Qu.: 6041233   1st Qu.: 6.0412  
##  Median : 8845293   Median : 8.8453  
##  Mean   : 9012886   Mean   : 9.0129  
##  3rd Qu.:11994827   3rd Qu.:11.9948  
##  Max.   :20262664   Max.   :20.2627
#Pivot Table To ID Top Fertilizer Applications by Country Sum 
sorghum.pivot<-dcast(data, COUNTRY ~., value.var="sorghum_Fe", fun.aggregate=sum)
sorghum.pivot<-arrange(sorghum.pivot,.)
tail(sorghum.pivot, 3) #ID Top 3 Fertilizer Countries
##           COUNTRY        .
## 161        Mexico 10571.98
## 162 United States 28892.45
## 163         China 68137.94
sum(sorghum.pivot$.) #Total Fertilizer
## [1] 217171.2
#Filter Countries
china<-filter(data, COUNTRY == "China") 
china.Perc<-sum(china$sorghum_Fe)/sum(sorghum.pivot$.) 
print(china.Perc) # % of Global Fertilizer For Crop 
## [1] 0.3137523
china.color<-"#1665AF"
summary(china)
##    OBJECTID           COUNTRY            Fishnet_ID      Fishnet__1   
##  Length:1417        Length:1417        Min.   :  464   Min.   :  464  
##  Class :character   Class :character   1st Qu.:48702   1st Qu.:48702  
##  Mode  :character   Mode  :character   Median :57638   Median :57638  
##                                        Mean   :54982   Mean   :54982  
##                                        3rd Qu.:64219   3rd Qu.:64219  
##                                        Max.   :76018   Max.   :76018  
##    barley_Fer        cassava_Fe        groundnut_         maize_Fert      
##  Min.   : 0.3023   Min.   : 0.5672   Min.   :  0.5672   Min.   :  0.5672  
##  1st Qu.:33.6040   1st Qu.:36.2767   1st Qu.: 30.5924   1st Qu.: 52.4740  
##  Median :35.4694   Median :37.5641   Median : 40.0470   Median : 66.0109  
##  Mean   :36.9683   Mean   :39.1958   Mean   : 40.7972   Mean   : 61.0785  
##  3rd Qu.:44.4152   3rd Qu.:47.9056   3rd Qu.: 52.7165   3rd Qu.: 75.1105  
##  Max.   :92.0236   Max.   :62.6335   Max.   :126.9427   Max.   :227.1572  
##    millet_Fer         rapeseed_F         rice_FertM         rye_FertMe     
##  Min.   :  0.5672   Min.   :  0.5672   Min.   :  0.4882   Min.   : 0.5672  
##  1st Qu.:  8.0391   1st Qu.: 54.1873   1st Qu.: 26.4062   1st Qu.:31.9068  
##  Median : 33.5524   Median : 58.4773   Median : 62.7289   Median :35.4209  
##  Mean   : 30.2364   Mean   : 60.2659   Mean   : 56.8110   Mean   :36.6549  
##  3rd Qu.: 44.1182   3rd Qu.: 71.7545   3rd Qu.: 79.1286   3rd Qu.:42.6427  
##  Max.   :105.8587   Max.   :151.1022   Max.   :158.1863   Max.   :92.6647  
##    sorghum_Fe         soybean_Fe         sunflower_        wheat_Fert      
##  Min.   :  0.5672   Min.   :  0.5672   Min.   : 0.5672   Min.   :  0.5481  
##  1st Qu.: 39.7837   1st Qu.: 30.7152   1st Qu.:18.3559   1st Qu.: 52.8928  
##  Median : 51.7188   Median : 40.5146   Median :37.4452   Median : 59.8756  
##  Mean   : 48.0861   Mean   : 40.9564   Mean   :33.6039   Mean   : 61.0944  
##  3rd Qu.: 66.6392   3rd Qu.: 49.7050   3rd Qu.:47.3547   3rd Qu.: 71.2560  
##  Max.   :116.3278   Max.   :111.8416   Max.   :89.3264   Max.   :148.6335  
##    distances          rescale_ND    
##  Min.   : 9387744   Min.   : 9.388  
##  1st Qu.:11428630   1st Qu.:11.429  
##  Median :12695272   Median :12.695  
##  Mean   :12671129   Mean   :12.671  
##  3rd Qu.:13863765   3rd Qu.:13.864  
##  Max.   :15943829   Max.   :15.944
usa<-filter(data, COUNTRY == "United States")
usa.Perc<-sum(usa$sorghum_Fe)/sum(sorghum.pivot$.) 
print(usa.Perc) # % of Global Fertilizer For Crop 
## [1] 0.13304
usa.color<-"#8F4A33"

mexico<-filter(data, COUNTRY == "Mexico")
mexico.Perc<-sum(mexico$sorghum_Fe)/sum(sorghum.pivot$.) 
print(mexico.Perc) # % of Global Fertilizer For Crop 
## [1] 0.04868041
mexico.color<-"#96AEEA"

#Plot 
ggplot(data, aes(x=rescale_ND,y=log10(data$sorghum_Fe))) + 
      geom_point(alpha=0.3) +theme_justin + 
      geom_point(data=china, aes(x=rescale_ND, y=log10(china$sorghum_Fe)), color=china.color, alpha=0.1) + 
      geom_point(data=usa, aes(x=rescale_ND, y=log10(usa$sorghum_Fe)), color=usa.color, alpha=0.2) + 
      geom_point(data=mexico, aes(x=rescale_ND, y=log10(mexico$sorghum_Fe)), color=mexico.color, alpha=0.1) + 
      geom_smooth(method="lm", color="red") +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Fertilizer Input")

###  Model - Sorghum Fertilizer Near
library(nlme)
data$logSorghum<-log10(data$sorghum_Fe)
# regular OLS no variance structure
sorghum.ols <- gls(logSorghum ~ rescale_ND, data = data)
# varFixed (variance changes linearly with X)
sorghum.fixed <- update(sorghum.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
sorghum.power <- update(sorghum.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
#sorghum.exp <- update(sorghum.ols,.~., weights = varExp(form = ~ rescale_ND)) # does not converge
# varConstPower (constant plus a power function of X (useful if X includes 0))
sorghum.ConstPower <- update(sorghum.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(sorghum.ols, sorghum.fixed, sorghum.power, sorghum.ConstPower) #sorghum.exp)
##                    df      AIC
## sorghum.ols         3 29910.81
## sorghum.fixed       3 30388.32
## sorghum.power       4 29412.24
## sorghum.ConstPower  5 29259.77
summary(sorghum.ConstPower)
## Generalized least squares fit by REML
##   Model: logSorghum ~ rescale_ND 
##   Data: data 
##        AIC      BIC    logLik
##   29259.77 29298.12 -14624.89
## 
## Variance function:
##  Structure: Constant plus power of variance covariate
##  Formula: ~rescale_ND 
##  Parameter estimates:
##      const      power 
## 696.211604   2.200842 
## 
## Coefficients:
##                 Value   Std.Error  t-value p-value
## (Intercept) 0.6204376 0.011306008 54.87681       0
## rescale_ND  0.0122760 0.001283534  9.56424       0
## 
##  Correlation: 
##            (Intr)
## rescale_ND -0.908
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -2.82904588 -0.84457116 -0.01915802  0.86712275  4.45485460 
## 
## Residual standard error: 0.0007211843 
## Degrees of freedom: 15830 total; 15828 residual

Soybean

soybean<- readOGR(dsn= "/Users/justinstewart/Dropbox/Collaborations/TobyKiers/Crop_Productivity_GCO/Analysis/GIS/Near_Fertilizer/", layer="Soy_Millet_Centroid") #Import GCO Centroid
## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/justinstewart/Dropbox/Collaborations/TobyKiers/Crop_Productivity_GCO/Analysis/GIS/Near_Fertilizer", layer: "Soy_Millet_Centroid"
## with 1 features
## It has 6 fields
## Integer64 fields read as strings:  OBJECTID ORIG_FID
soybean.c<-coordinates(soybean) #Extract Lat-Long
distances <- spDistsN1(xy,soybean.c, longlat = FALSE) #Calculate Near Distances
data<-cbind(data.frame(metadata),data.frame(distances)) #Merge into a DF
data$rescale_ND<-data$distances/(1000*1000) #Rescale distances 
data<-data %>% filter(soybean_Fe > 0, na.rm=TRUE) #select Fertilizer > 0 

#Pivot Table To ID Top Fertilizer Applications by Country Sum 
soybean.pivot<-dcast(data, COUNTRY ~., value.var="soybean_Fe", fun.aggregate=sum)
soybean.pivot<-arrange(soybean.pivot,.)
tail(soybean.pivot, 3) #ID Top 3 Fertilizer Countries
##           COUNTRY        .
## 161        Brazil 16426.22
## 162 United States 24062.62
## 163         China 58035.28
sum(soybean.pivot$.) #Total Fertilizer
## [1] 192474.3
#Filter Countries
china<-filter(data, COUNTRY == "China") 
china.Perc<-sum(china$soybean_Fe)/sum(soybean.pivot$.) 
print(china.Perc) # % of Global Fertilizer For Crop 
## [1] 0.3015222
china.color<-"#1665AF"
summary(china)
##    OBJECTID           COUNTRY            Fishnet_ID      Fishnet__1   
##  Length:1417        Length:1417        Min.   :  464   Min.   :  464  
##  Class :character   Class :character   1st Qu.:48702   1st Qu.:48702  
##  Mode  :character   Mode  :character   Median :57638   Median :57638  
##                                        Mean   :54982   Mean   :54982  
##                                        3rd Qu.:64219   3rd Qu.:64219  
##                                        Max.   :76018   Max.   :76018  
##    barley_Fer        cassava_Fe        groundnut_         maize_Fert      
##  Min.   : 0.3023   Min.   : 0.5672   Min.   :  0.5672   Min.   :  0.5672  
##  1st Qu.:33.6040   1st Qu.:36.2767   1st Qu.: 30.5924   1st Qu.: 52.4740  
##  Median :35.4694   Median :37.5641   Median : 40.0470   Median : 66.0109  
##  Mean   :36.9683   Mean   :39.1958   Mean   : 40.7972   Mean   : 61.0785  
##  3rd Qu.:44.4152   3rd Qu.:47.9056   3rd Qu.: 52.7165   3rd Qu.: 75.1105  
##  Max.   :92.0236   Max.   :62.6335   Max.   :126.9427   Max.   :227.1572  
##    millet_Fer         rapeseed_F         rice_FertM         rye_FertMe     
##  Min.   :  0.5672   Min.   :  0.5672   Min.   :  0.4882   Min.   : 0.5672  
##  1st Qu.:  8.0391   1st Qu.: 54.1873   1st Qu.: 26.4062   1st Qu.:31.9068  
##  Median : 33.5524   Median : 58.4773   Median : 62.7289   Median :35.4209  
##  Mean   : 30.2364   Mean   : 60.2659   Mean   : 56.8110   Mean   :36.6549  
##  3rd Qu.: 44.1182   3rd Qu.: 71.7545   3rd Qu.: 79.1286   3rd Qu.:42.6427  
##  Max.   :105.8587   Max.   :151.1022   Max.   :158.1863   Max.   :92.6647  
##    sorghum_Fe         soybean_Fe         sunflower_        wheat_Fert      
##  Min.   :  0.5672   Min.   :  0.5672   Min.   : 0.5672   Min.   :  0.5481  
##  1st Qu.: 39.7837   1st Qu.: 30.7152   1st Qu.:18.3559   1st Qu.: 52.8928  
##  Median : 51.7188   Median : 40.5146   Median :37.4452   Median : 59.8756  
##  Mean   : 48.0861   Mean   : 40.9564   Mean   :33.6039   Mean   : 61.0944  
##  3rd Qu.: 66.6392   3rd Qu.: 49.7050   3rd Qu.:47.3547   3rd Qu.: 71.2560  
##  Max.   :116.3278   Max.   :111.8416   Max.   :89.3264   Max.   :148.6335  
##    distances          rescale_ND    
##  Min.   : 9387671   Min.   : 9.388  
##  1st Qu.:11428558   1st Qu.:11.429  
##  Median :12695199   Median :12.695  
##  Mean   :12671056   Mean   :12.671  
##  3rd Qu.:13863692   3rd Qu.:13.864  
##  Max.   :15943756   Max.   :15.944
usa<-filter(data, COUNTRY == "United States")
usa.Perc<-sum(usa$soybean_Fe)/sum(soybean.pivot$.) 
print(usa.Perc) # % of Global Fertilizer For Crop 
## [1] 0.1250173
usa.color<-"#8F4A33"

brazil<-filter(data, COUNTRY == "Brazil")
brazil.Perc<-sum(brazil$soybean_Fe)/sum(cassava.pivot$.) 
print(brazil.Perc) # % of Global Fertilizer For Crop 
## [1] 0.1060901
brazil.color<-"#E28D15"

#Plot 
ggplot(data, aes(x=rescale_ND,y=log10(data$soybean_Fe))) + 
      geom_point(alpha=0.3) +theme_justin + 
      geom_point(data=china, aes(x=rescale_ND, y=log10(china$soybean_Fe)), color=china.color, alpha=0.1) + 
      geom_point(data=usa, aes(x=rescale_ND, y=log10(usa$soybean_Fe)), color=usa.color, alpha=0.2) + 
      geom_point(data=brazil, aes(x=rescale_ND, y=log10(brazil$soybean_Fe)), color=brazil.color, alpha=0.3) + 
      geom_smooth(method="lm", color="red") +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Fertilizer Input")

###  Model - Soybean Fertilizer Near
library(nlme)
data$logSoybean<-log10(data$soybean_Fe)
# regular OLS no variance structure
soybean.ols <- gls(logSoybean ~ rescale_ND, data = data)
# varFixed (variance changes linearly with X)
soybean.fixed <- update(soybean.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
soybean.power <- update(soybean.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
soybean.exp <- update(soybean.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
soybean.ConstPower <- update(soybean.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(soybean.ols, soybean.fixed, soybean.power, soybean.ConstPower, soybean.exp)
##                    df      AIC
## soybean.ols         3 28705.23
## soybean.fixed       3 27951.06
## soybean.power       4 27776.14
## soybean.ConstPower  5 27778.14
## soybean.exp         4 27937.14
summary(soybean.exp)
## Generalized least squares fit by REML
##   Model: logSoybean ~ rescale_ND 
##   Data: data 
##        AIC      BIC    logLik
##   27937.14 27967.82 -13964.57
## 
## Variance function:
##  Structure: Exponential of variance covariate
##  Formula: ~rescale_ND 
##  Parameter estimates:
##      expon 
## 0.04240481 
## 
## Coefficients:
##                 Value   Std.Error t-value p-value
## (Intercept) 0.4743330 0.010283820 46.1242       0
## rescale_ND  0.0258654 0.001192401 21.6919       0
## 
##  Correlation: 
##            (Intr)
## rescale_ND -0.898
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -2.8424926 -0.7432165 -0.1372582  0.8539977  2.8222918 
## 
## Residual standard error: 0.3986962 
## Degrees of freedom: 15830 total; 15828 residual

Wheat

wheat<- readOGR(dsn= "/Users/justinstewart/Dropbox/Collaborations/TobyKiers/Crop_Productivity_GCO/Analysis/GIS/Near_Fertilizer/", layer="Barley_Wheat_Centroid") #Import GCO Centroid
## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/justinstewart/Dropbox/Collaborations/TobyKiers/Crop_Productivity_GCO/Analysis/GIS/Near_Fertilizer", layer: "Barley_Wheat_Centroid"
## with 1 features
## It has 2 fields
## Integer64 fields read as strings:  OBJECTID ORIG_FID
wheat.c<-coordinates(wheat) #Extract Lat-Long
distances <- spDistsN1(xy,wheat.c, longlat = FALSE) #Calculate Near Distances
data<-cbind(data.frame(metadata),data.frame(distances)) #Merge into a DF
data$rescale_ND<-data$distances/(1000*1000) #Rescale distances 
data<-data %>% filter(wheat_Fert > 0, na.rm=TRUE) #select Fertilizer > 0 

#Pivot Table To ID Top Fertilizer Applications by Country Sum 
wheat.pivot<-dcast(data, COUNTRY ~., value.var="wheat_Fert", fun.aggregate=sum)
wheat.pivot<-arrange(wheat.pivot,.)
tail(wheat.pivot, 3) #ID Top 3 Fertilizer Countries
##           COUNTRY        .
## 161         India 19352.26
## 162 United States 42375.82
## 163         China 86570.81
sum(wheat.pivot$.) #Total Fertilizer
## [1] 389705
#Filter Countries
china<-filter(data, COUNTRY == "China") 
china.Perc<-sum(china$wheat_Fert)/sum(wheat.pivot$.) 
print(china.Perc) # % of Global Fertilizer For Crop 
## [1] 0.2221445
china.color<-"#1665AF"

usa<-filter(data, COUNTRY == "United States")
usa.Perc<-sum(usa$wheat_Fert)/sum(wheat.pivot$.) 
print(usa.Perc) # % of Global Fertilizer For Crop 
## [1] 0.1087382
usa.color<-"#8F4A33"

india<-filter(data, COUNTRY == "India")
india.Perc<-sum(canada$wheat_Fert)/sum(wheat.pivot$.) 
print(india.Perc) # % of Global Fertilizer For Crop 
## [1] 0.03751558
india.color<-"#078D40"

#Plot 
ggplot(data, aes(x=rescale_ND,y=log10(data$wheat_Fert))) + 
      geom_point(alpha=0.3) +theme_justin + 
      geom_point(data=china, aes(x=rescale_ND, y=log10(china$wheat_Fert)), color=china.color, alpha=0.1) + 
      geom_point(data=usa, aes(x=rescale_ND, y=log10(usa$wheat_Fert)), color=usa.color, alpha=0.2) + 
      geom_point(data=india, aes(x=rescale_ND, y=log10(india$wheat_Fert)), color=india.color, alpha=0.3) + 
      geom_smooth(method="lm", color="red") +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Fertilizer Input")

###  Model - Wheat Fertilizer Near
library(nlme)
data$logWheat<-log10(data$wheat_Fe)
# regular OLS no variance structure
wheat.ols <- gls(logWheat ~ rescale_ND, data = data)
# varFixed (variance changes linearly with X)
wheat.fixed <- update(wheat.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
wheat.power <- update(wheat.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
wheat.exp <- update(wheat.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
wheat.ConstPower <- update(wheat.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(wheat.ols, wheat.fixed, wheat.power, wheat.ConstPower, wheat.exp)
##                  df      AIC
## wheat.ols         3 30526.38
## wheat.fixed       3 31649.24
## wheat.power       4 30472.62
## wheat.ConstPower  5 30474.62
## wheat.exp         4 30474.90
summary(wheat.power)
## Generalized least squares fit by REML
##   Model: logWheat ~ rescale_ND 
##   Data: data 
##        AIC     BIC    logLik
##   30472.62 30503.3 -15232.31
## 
## Variance function:
##  Structure: Power of variance covariate
##  Formula: ~rescale_ND 
##  Parameter estimates:
##      power 
## 0.08836862 
## 
## Coefficients:
##                 Value   Std.Error  t-value p-value
## (Intercept) 0.6647638 0.012037325 55.22521       0
## rescale_ND  0.0398286 0.001266645 31.44416       0
## 
##  Correlation: 
##            (Intr)
## rescale_ND -0.909
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -3.0217546 -0.8158082  0.1384038  0.8436922  2.9742344 
## 
## Residual standard error: 0.527033 
## Degrees of freedom: 15830 total; 15828 residual

Maize

maize<- readOGR(dsn= "/Users/justinstewart/Dropbox/Collaborations/TobyKiers/Crop_Productivity_GCO/Analysis/GIS/Near_Fertilizer/", layer="Maize") #Import GCO Centroid
## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/justinstewart/Dropbox/Collaborations/TobyKiers/Crop_Productivity_GCO/Analysis/GIS/Near_Fertilizer", layer: "Maize"
## with 1 features
## It has 2 fields
## Integer64 fields read as strings:  OBJECTID ORIG_FID
maize.c<-coordinates(maize) #Extract Lat-Long
distances <- spDistsN1(xy,maize.c, longlat = FALSE) #Calculate Near Distances
data<-cbind(data.frame(metadata),data.frame(distances)) #Merge into a DF
data$rescale_ND<-data$distances/(1000*1000) #Rescale distances 
data<-data %>% filter(maize_Fert > 0, na.rm=TRUE) #select Fertilizer > 0 

#Pivot Table To ID Top Fertilizer Applications by Country Sum 
maize.pivot<-dcast(data, COUNTRY ~., value.var="maize_Fert", fun.aggregate=sum)
maize.pivot<-arrange(maize.pivot,.)
tail(maize.pivot, 3) #ID Top 3 Fertilizer Countries
##           COUNTRY        .
## 161        Brazil 26921.38
## 162 United States 71878.71
## 163         China 86548.20
sum(maize.pivot$.) #Total Fertilizer
## [1] 461308.8
#Filter Countries
china<-filter(data, COUNTRY == "China") 
china.Perc<-sum(china$maize_Fert)/sum(maize.pivot$.) 
print(china.Perc) # % of Global Fertilizer For Crop 
## [1] 0.1876144
china.color<-"#1665AF"

usa<-filter(data, COUNTRY == "United States")
usa.Perc<-sum(usa$maize_Fert)/sum(maize.pivot$.) 
print(usa.Perc) # % of Global Fertilizer For Crop 
## [1] 0.1558147
usa.color<-"#8F4A33"

brazil<-filter(data, COUNTRY == "Brazil")
brazil.Perc<-sum(brazil$maize_Fert)/sum(maize.pivot$.) 
print(brazil.Perc) # % of Global Fertilizer For Crop 
## [1] 0.0583587
brazil.color<-"#E28D15"

#Plot 
ggplot(data, aes(x=rescale_ND,y=log10(data$maize_Fert))) + 
      geom_point(alpha=0.3) +theme_justin + 
      geom_point(data=china, aes(x=rescale_ND, y=log10(china$maize_Fert)), color=china.color, alpha=0.1) + 
      geom_point(data=usa, aes(x=rescale_ND, y=log10(usa$maize_Fert)), color=usa.color, alpha=0.2) + 
      geom_point(data=brazil, aes(x=rescale_ND, y=log10(brazil$maize_Fert)), color=brazil.color, alpha=0.3) + 
      geom_smooth(method="lm", color="red") +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Fertilizer Input")

###  Model - Maize Fertilizer Near
library(nlme)
data$logMaize<-log10(data$maize_Fert)
# regular OLS no variance structure
maize.ols <- gls(logMaize ~ rescale_ND, data = data)
# varFixed (variance changes linearly with X)
maize.fixed <- update(maize.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
maize.power <- update(maize.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
maize.exp <- update(maize.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
maize.ConstPower <- update(maize.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(maize.ols, maize.fixed, maize.power, maize.ConstPower, maize.exp)
##                  df      AIC
## maize.ols         3 31362.07
## maize.fixed       3 32683.88
## maize.power       4 31207.87
## maize.ConstPower  5 31099.31
## maize.exp         4 31149.45
summary(maize.ConstPower)
## Generalized least squares fit by REML
##   Model: logMaize ~ rescale_ND 
##   Data: data 
##        AIC      BIC    logLik
##   31099.31 31137.65 -15544.65
## 
## Variance function:
##  Structure: Constant plus power of variance covariate
##  Formula: ~rescale_ND 
##  Parameter estimates:
##       const       power 
## 8370.495676    2.896854 
## 
## Coefficients:
##                 Value   Std.Error  t-value p-value
## (Intercept) 0.8558426 0.012497054 68.48355       0
## rescale_ND  0.0272339 0.001361917 19.99677       0
## 
##  Correlation: 
##            (Intr)
## rescale_ND -0.913
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -3.1721489 -0.6182477  0.1404735  0.8187922  3.2967145 
## 
## Residual standard error: 7.005263e-05 
## Degrees of freedom: 15830 total; 15828 residual

Millet

millet<- readOGR(dsn= "/Users/justinstewart/Dropbox/Collaborations/TobyKiers/Crop_Productivity_GCO/Analysis/GIS/Near_Fertilizer/", layer="Soy_Millet_Centroid") #Import GCO Centroid
## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/justinstewart/Dropbox/Collaborations/TobyKiers/Crop_Productivity_GCO/Analysis/GIS/Near_Fertilizer", layer: "Soy_Millet_Centroid"
## with 1 features
## It has 6 fields
## Integer64 fields read as strings:  OBJECTID ORIG_FID
millet.c<-coordinates(millet) #Extract Lat-Long
distances <- spDistsN1(xy,millet.c, longlat = FALSE) #Calculate Near Distances
data<-cbind(data.frame(metadata),data.frame(distances)) #Merge into a DF
data$rescale_ND<-data$distances/(1000*1000) #Rescale distances 
data<-data %>% filter(millet_Fer > 0, na.rm=TRUE) #select Fertilizer > 0 

#Pivot Table To ID Top Fertilizer Applications by Country Sum 
millet.pivot<-dcast(data, COUNTRY ~., value.var="millet_Fer", fun.aggregate=sum)
millet.pivot<-arrange(millet.pivot,.)
tail(millet.pivot, 3) #ID Top 3 Fertilizer Countries
##           COUNTRY         .
## 161         India  8254.568
## 162 United States  8801.661
## 163         China 42844.987
sum(millet.pivot$.) #Total Fertilizer
## [1] 140288.2
#Filter Countries
china<-filter(data, COUNTRY == "China") 
china.Perc<-sum(china$millet_Fer)/sum(millet.pivot$.) 
print(china.Perc) # % of Global Fertilizer For Crop 
## [1] 0.3054069
china.color<-"#1665AF"

usa<-filter(data, COUNTRY == "United States")
usa.Perc<-sum(usa$millet_Fer)/sum(millet.pivot$.) 
print(usa.Perc) # % of Global Fertilizer For Crop 
## [1] 0.06273986
usa.color<-"#8F4A33"

india<-filter(data, COUNTRY == "India")
india.Perc<-sum(india$millet_Fer)/sum(millet.pivot$.) 
print(india.Perc) # % of Global Fertilizer For Crop 
## [1] 0.05884007
india.color<-"#078D40"

#Plot 
ggplot(data, aes(x=rescale_ND,y=log10(data$millet_Fer))) + 
      geom_point(alpha=0.3) +theme_justin + 
      geom_point(data=china, aes(x=rescale_ND, y=log10(china$millet_Fer)), color=china.color, alpha=0.1) + 
      geom_point(data=usa, aes(x=rescale_ND, y=log10(usa$millet_Fer)), color=usa.color, alpha=0.2) + 
      geom_point(data=india, aes(x=rescale_ND, y=log10(india$millet_Fer)), color=india.color, alpha=0.3) + 
      geom_smooth(method="lm", color="red") +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Fertilizer Input")

###  Model - Millet Fertilizer Near
library(nlme)
data$logMillet<-log10(data$millet_Fer)
# regular OLS no variance structure
millet.ols <- gls(logMillet ~ rescale_ND, data = data)
# varFixed (variance changes linearly with X)
millet.fixed <- update(millet.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
millet.power <- update(millet.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
millet.exp <- update(millet.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
millet.ConstPower <- update(millet.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(millet.ols, millet.fixed, millet.power, millet.ConstPower, millet.exp)
##                   df      AIC
## millet.ols         3 26420.45
## millet.fixed       3 27286.77
## millet.power       4 26017.22
## millet.ConstPower  5 25576.64
## millet.exp         4 25803.51
summary(millet.ConstPower)
## Generalized least squares fit by REML
##   Model: logMillet ~ rescale_ND 
##   Data: data 
##        AIC      BIC    logLik
##   25576.64 25614.98 -12783.32
## 
## Variance function:
##  Structure: Constant plus power of variance covariate
##  Formula: ~rescale_ND 
##  Parameter estimates:
##        const        power 
## 1.043281e+05 4.031497e+00 
## 
## Coefficients:
##                  Value   Std.Error t-value p-value
## (Intercept)  0.5954364 0.010446564 56.9983  0.0000
## rescale_ND  -0.0003623 0.001180226 -0.3070  0.7588
## 
##  Correlation: 
##            (Intr)
## rescale_ND -0.915
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -3.0304616 -0.7617620 -0.2025050  0.7175356  3.4173822 
## 
## Residual standard error: 4.568136e-06 
## Degrees of freedom: 15830 total; 15828 residual

Rapeseed

rapeseed<- readOGR(dsn= "/Users/justinstewart/Dropbox/Collaborations/TobyKiers/Crop_Productivity_GCO/Analysis/GIS/Near_Fertilizer/", layer="Rapeseed_Centroid") #Import GCO Centroid
## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/justinstewart/Dropbox/Collaborations/TobyKiers/Crop_Productivity_GCO/Analysis/GIS/Near_Fertilizer", layer: "Rapeseed_Centroid"
## with 1 features
## It has 2 fields
## Integer64 fields read as strings:  OBJECTID ORIG_FID
rapeseed.c<-coordinates(rapeseed) #Extract Lat-Long
distances <- spDistsN1(xy,rapeseed.c, longlat = FALSE) #Calculate Near Distances
data<-cbind(data.frame(metadata),data.frame(distances)) #Merge into a DF
data$rescale_ND<-data$distances/(1000*1000) #Rescale distances 
data<-data %>% filter(rapeseed_F > 0, na.rm=TRUE) #select Fertilizer > 0 

#Pivot Table To ID Top Fertilizer Applications by Country Sum 
rapeseed.pivot<-dcast(data, COUNTRY ~., value.var="rapeseed_F", fun.aggregate=sum)
rapeseed.pivot<-arrange(rapeseed.pivot,.)
tail(rapeseed.pivot, 3) #ID Top 3 Fertilizer Countries
##           COUNTRY        .
## 161 United States 23782.47
## 162        Brazil 27471.87
## 163         China 85396.80
sum(rapeseed.pivot$.) #Total Fertilizer
## [1] 331872.4
#Filter Countries
china<-filter(data, COUNTRY == "China") 
china.Perc<-sum(china$rapeseed_F)/sum(rapeseed.pivot$.) 
print(china.Perc) # % of Global Fertilizer For Crop 
## [1] 0.2573182
china.color<-"#1665AF"

brazil<-filter(data, COUNTRY == "Brazil")
brazil.Perc<-sum(brazil$rapeseed_F)/sum(rapeseed.pivot$.) 
print(brazil.Perc) # % of Global Fertilizer For Crop 
## [1] 0.08277842
brazil.color<-"#E28D15"

usa<-filter(data, COUNTRY == "United States")
usa.Perc<-sum(usa$rapeseed_F)/sum(rapeseed.pivot$.) 
print(usa.Perc) # % of Global Fertilizer For Crop 
## [1] 0.07166151
usa.color<-"#8F4A33"

#Plot 
ggplot(data, aes(x=rescale_ND,y=log10(data$rapeseed_F))) + 
      geom_point(alpha=0.3) +theme_justin + 
      geom_point(data=china, aes(x=rescale_ND, y=log10(china$rapeseed_F)), color=china.color, alpha=0.1) + 
      geom_point(data=usa, aes(x=rescale_ND, y=log10(usa$rapeseed_F)), color=usa.color, alpha=0.2) + 
      geom_point(data=brazil, aes(x=rescale_ND, y=log10(brazil$rapeseed_F)), color=brazil.color, alpha=0.3) + 
      geom_smooth(method="lm", color="red") +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Fertilizer Input")

###  Model - Rapeseed Fertilizer Near
library(nlme)
data$logRapeseed<-log10(data$rapeseed_F)
# regular OLS no variance structure
rapeseed.ols <- gls(logRapeseed ~ rescale_ND, data = data)
# varFixed (variance changes linearly with X)
rapeseed.fixed <- update(rapeseed.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
rapeseed.power <- update(rapeseed.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
rapeseed.exp <- update(rapeseed.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
rapeseed.ConstPower <- update(rapeseed.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(rapeseed.ols, rapeseed.fixed, rapeseed.power, rapeseed.ConstPower, rapeseed.exp)
##                     df      AIC
## rapeseed.ols         3 31695.21
## rapeseed.fixed       3 35106.56
## rapeseed.power       4 31657.72
## rapeseed.ConstPower  5 31415.93
## rapeseed.exp         4 31496.68
summary(rapeseed.ConstPower)
## Generalized least squares fit by REML
##   Model: logRapeseed ~ rescale_ND 
##   Data: data 
##        AIC      BIC    logLik
##   31415.93 31454.27 -15702.96
## 
## Variance function:
##  Structure: Constant plus power of variance covariate
##  Formula: ~rescale_ND 
##  Parameter estimates:
##       const       power 
## 2998.962463    2.522422 
## 
## Coefficients:
##                 Value   Std.Error   t-value p-value
## (Intercept)  1.008511 0.012178356  82.81171       0
## rescale_ND  -0.014636 0.001284464 -11.39466       0
## 
##  Correlation: 
##            (Intr)
## rescale_ND -0.907
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -2.6206202 -0.8438772 -0.2284213  0.9625367  2.0841019 
## 
## Residual standard error: 0.0001939111 
## Degrees of freedom: 15830 total; 15828 residual

Rice

rice<- readOGR(dsn= "/Users/justinstewart/Dropbox/Collaborations/TobyKiers/Crop_Productivity_GCO/Analysis/GIS/Near_Fertilizer/", layer="Rice_Centroid") #Import GCO Centroid
## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/justinstewart/Dropbox/Collaborations/TobyKiers/Crop_Productivity_GCO/Analysis/GIS/Near_Fertilizer", layer: "Rice_Centroid"
## with 1 features
## It has 2 fields
## Integer64 fields read as strings:  OBJECTID ORIG_FID
rice.c<-coordinates(rice) #Extract Lat-Long
distances <- spDistsN1(xy,rice.c, longlat = FALSE) #Calculate Near Distances
data<-cbind(data.frame(metadata),data.frame(distances)) #Merge into a DF
data$rescale_ND<-data$distances/(1000*1000) #Rescale distances 
data<-data %>% filter(rice_FertM > 0, na.rm=TRUE) #select Fertilizer > 0 

#Pivot Table To ID Top Fertilizer Applications by Country Sum 
rice.pivot<-dcast(data, COUNTRY ~., value.var="rice_FertM", fun.aggregate=sum)
rice.pivot<-arrange(rice.pivot,.)
tail(rice.pivot, 3) #ID Top 3 Fertilizer Countries
##     COUNTRY        .
## 161  Brazil 16411.96
## 162   India 19194.99
## 163   China 80501.18
sum(rice.pivot$.) #Total Fertilizer
## [1] 261110.5
#Filter Countries
china<-filter(data, COUNTRY == "China") 
china.Perc<-sum(china$rice_FertM)/sum(rice.pivot$.) 
print(china.Perc) # % of Global Fertilizer For Crop 
## [1] 0.3083031
china.color<-"#1665AF"

india<-filter(data, COUNTRY == "India")
india.Perc<-sum(india$rice_FertM)/sum(rice.pivot$.) 
print(india.Perc) # % of Global Fertilizer For Crop 
## [1] 0.07351288
india.color<-"#078D40"

brazil<-filter(data, COUNTRY == "Brazil")
brazil.Perc<-sum(brazil$rice_FertM)/sum(rice.pivot$.) 
print(brazil.Perc) # % of Global Fertilizer For Crop 
## [1] 0.06285444
brazil.color<-"#E28D15"

#Plot 
ggplot(data, aes(x=rescale_ND,y=log10(data$rice_Fert))) + 
      geom_point(alpha=0.3) +theme_justin + 
      geom_point(data=china, aes(x=rescale_ND, y=log10(china$rice_Fert)), color=china.color, alpha=0.1) + 
      geom_point(data=india, aes(x=rescale_ND, y=log10(india$rice_Fert)), color=india.color, alpha=0.2) + 
      geom_point(data=brazil, aes(x=rescale_ND, y=log10(brazil$rice_Fert)), color=brazil.color, alpha=0.3) + 
      geom_smooth(method="lm", color="red") +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Fertilizer Input")

###  Model - Rice Fertilizer Near
library(nlme)
data$logRice<-log10(data$rice_Fert)
# regular OLS no variance structure
rice.ols <- gls(logRice ~ rescale_ND, data = data)
# varFixed (variance changes linearly with X)
rice.fixed <- update(rice.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
rice.power <- update(rice.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
rice.exp <- update(rice.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
rice.ConstPower <- update(rice.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(rice.ols, rice.fixed, rice.power, rice.ConstPower, rice.exp)
##                 df      AIC
## rice.ols         3 31997.48
## rice.fixed       3 31055.43
## rice.power       4 30845.14
## rice.ConstPower  5 30701.07
## rice.exp         4 30705.02
summary(rice.ConstPower)
## Generalized least squares fit by REML
##   Model: logRice ~ rescale_ND 
##   Data: data 
##        AIC      BIC    logLik
##   30701.07 30739.42 -15345.54
## 
## Variance function:
##  Structure: Constant plus power of variance covariate
##  Formula: ~rescale_ND 
##  Parameter estimates:
##    const    power 
## 80.82712  1.63595 
## 
## Coefficients:
##                 Value   Std.Error  t-value p-value
## (Intercept) 0.6502007 0.010969286 59.27466       0
## rescale_ND  0.0126999 0.001317242  9.64124       0
## 
##  Correlation: 
##            (Intr)
## rescale_ND -0.897
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -2.8569811 -0.8453237 -0.1289194  0.9090870  2.6058951 
## 
## Residual standard error: 0.005390671 
## Degrees of freedom: 15830 total; 15828 residual

Rye

rye<- readOGR(dsn= "/Users/justinstewart/Dropbox/Collaborations/TobyKiers/Crop_Productivity_GCO/Analysis/GIS/Near_Fertilizer/", layer="Rye_Centroid") #Import GCO Centroid
## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/justinstewart/Dropbox/Collaborations/TobyKiers/Crop_Productivity_GCO/Analysis/GIS/Near_Fertilizer", layer: "Rye_Centroid"
## with 1 features
## It has 2 fields
## Integer64 fields read as strings:  OBJECTID ORIG_FID
rye.c<-coordinates(rye) #Extract Lat-Long
distances <- spDistsN1(xy,rye.c, longlat = FALSE) #Calculate Near Distances
data<-cbind(data.frame(metadata),data.frame(distances)) #Merge into a DF
data$rescale_ND<-data$distances/(1000*1000) #Rescale distances 
data<-data %>% filter(rye_FertMe > 0, na.rm=TRUE) #select Fertilizer > 0 

#Pivot Table To ID Top Fertilizer Applications by Country Sum 
rye.pivot<-dcast(data, COUNTRY ~., value.var="rye_FertMe", fun.aggregate=sum)
rye.pivot<-arrange(rye.pivot,.)
tail(rye.pivot, 3) #ID Top 3 Fertilizer Countries
##                COUNTRY        .
## 161      United States 12028.64
## 162 Russian Federation 17563.30
## 163              China 51940.00
sum(rye.pivot$.) #Total Fertilizer
## [1] 199071.4
#Filter Countries
china<-filter(data, COUNTRY == "China") 
china.Perc<-sum(china$rye_Fert)/sum(rye.pivot$.) 
print(china.Perc) # % of Global Fertilizer For Crop 
## [1] 0.2609114
china.color<-"#1665AF"

usa<-filter(data, COUNTRY == "United States")
usa.Perc<-sum(usa$rye_Fert)/sum(rye.pivot$.) 
print(usa.Perc) # % of Global Fertilizer For Crop 
## [1] 0.06042376
usa.color<-"#8F4A33"

russia<-filter(data, COUNTRY == "Russian Federation")
russia.Perc<-sum(russia$rye_FertMe)/sum(rye.pivot$.) 
print(russia.Perc) # % of Global Fertilizer For Crop 
## [1] 0.08822614
russia.color<-"#D3867A"

#Plot 
ggplot(data, aes(x=rescale_ND,y=log10(data$rye_FertMe))) + 
      geom_point(alpha=0.3) +theme_justin + 
      geom_point(data=china, aes(x=rescale_ND, y=log10(china$rye_FertMe)), color=china.color, alpha=0.1) + 
      geom_point(data=usa, aes(x=rescale_ND, y=log10(usa$rye_FertMe)), color=usa.color, alpha=0.2) + 
      geom_point(data=russia, aes(x=rescale_ND, y=log10(russia$rye_FertMe)), color=russia.color, alpha=0.3) + 
      geom_smooth(method="lm", color="red") +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Fertilizer Input")

###  Model - Rye Fertilizer Near
library(nlme)
data$logRye<-log10(data$rye_FertMe)
# regular OLS no variance structure
rye.ols <- gls(logRye ~ rescale_ND, data = data)
# varFixed (variance changes linearly with X)
rye.fixed <- update(rye.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
rye.power <- update(rye.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
rye.exp <- update(rye.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
rye.ConstPower <- update(rye.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(rye.ols, rye.fixed, rye.power, rye.ConstPower, rye.exp)
##                df      AIC
## rye.ols         3 26797.09
## rye.fixed       3 32323.56
## rye.power       4 26789.40
## rye.ConstPower  5 26791.40
## rye.exp         4 26723.11
summary(rye.exp)
## Generalized least squares fit by REML
##   Model: logRye ~ rescale_ND 
##   Data: data 
##        AIC      BIC    logLik
##   26723.11 26753.79 -13357.56
## 
## Variance function:
##  Structure: Exponential of variance covariate
##  Formula: ~rescale_ND 
##  Parameter estimates:
##        expon 
## -0.009075993 
## 
## Coefficients:
##                  Value   Std.Error   t-value p-value
## (Intercept)  0.8319736 0.008414108  98.87841       0
## rescale_ND  -0.0090626 0.000767067 -11.81465       0
## 
##  Correlation: 
##            (Intr)
## rescale_ND -0.848
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -2.6996057 -0.8100850 -0.1754591  0.9082522  2.5407947 
## 
## Residual standard error: 0.6084994 
## Degrees of freedom: 15830 total; 15828 residual

Sunflower

sunflower<- readOGR(dsn= "/Users/justinstewart/Dropbox/Collaborations/TobyKiers/Crop_Productivity_GCO/Analysis/GIS/Near_Fertilizer/", layer="Sunflower_Centroid") #Import GCO Centroid
## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/justinstewart/Dropbox/Collaborations/TobyKiers/Crop_Productivity_GCO/Analysis/GIS/Near_Fertilizer", layer: "Sunflower_Centroid"
## with 1 features
## It has 2 fields
## Integer64 fields read as strings:  OBJECTID ORIG_FID
sunflower.c<-coordinates(sunflower) #Extract Lat-Long
distances <- spDistsN1(xy,sunflower.c, longlat = FALSE) #Calculate Near Distances
data<-cbind(data.frame(metadata),data.frame(distances)) #Merge into a DF
data$rescale_ND<-data$distances/(1000*1000) #Rescale distances 
data<-data %>% filter(sunflower_ > 0, na.rm=TRUE) #select Fertilizer > 0 

#Pivot Table To ID Top Fertilizer Applications by Country Sum 
sunflower.pivot<-dcast(data, COUNTRY ~., value.var="sunflower_", fun.aggregate=sum)
sunflower.pivot<-arrange(sunflower.pivot,.)
tail(sunflower.pivot, 3) #ID Top 3 Fertilizer Countries
##     COUNTRY        .
## 161   India 11277.46
## 162  France 13535.91
## 163   China 47616.79
sum(sunflower.pivot$.) #Total Fertilizer
## [1] 190902.9
#Filter Countries
china<-filter(data, COUNTRY == "China") 
china.Perc<-sum(china$sunflower_)/sum(sunflower.pivot$.) 
print(china.Perc) # % of Global Fertilizer For Crop 
## [1] 0.2494294
china.color<-"#1665AF"

france<-filter(data, COUNTRY == "France")
france.Perc<-sum(usa$sunflower_)/sum(sunflower.pivot$.) 
print(france.Perc) # % of Global Fertilizer For Crop 
## [1] 0.05220243
france.color<-"#CE465B"

india<-filter(data, COUNTRY == "India")
india.Perc<-sum(india$sunflower_)/sum(rice.pivot$.) 
print(india.Perc) # % of Global Fertilizer For Crop 
## [1] 0.04319038
india.color<-"#078D40"

#Plot 
ggplot(data, aes(x=rescale_ND,y=log10(data$sunflower_))) + 
      geom_point(alpha=0.3) +theme_justin + 
      geom_point(data=china, aes(x=rescale_ND, y=log10(china$sunflower_)), color=china.color, alpha=0.1) + 
      geom_point(data=france, aes(x=rescale_ND, y=log10(france$sunflower_)), color=france.color, alpha=0.2) + 
      geom_point(data=india, aes(x=rescale_ND, y=log10(india$sunflower_)), color=india.color, alpha=0.3) + 
      geom_smooth(method="lm", color="red") +xlab("Distance From GCO (1000's km)") +ylab ("Log10 Fertilizer Input")

###  Model - Sunflower Fertilizer Near
library(nlme)
data$logSunflower<-log10(data$sunflower_)
# regular OLS no variance structure
sunflower.ols <- gls(logSunflower ~ rescale_ND, data = data)
# varFixed (variance changes linearly with X)
sunflower.fixed <- update(sunflower.ols, .~., weights = varFixed(~rescale_ND))
# varPower (variance changes as a power function with X)
sunflower.power <- update(sunflower.ols, . ~ ., weights = varPower(form = ~rescale_ND))
# varExp (variance changes as an exponential function of x)
sunflower.exp <- update(sunflower.ols, . ~., weights = varExp(form = ~ rescale_ND)) # error
# varConstPower (constant plus a power function of X (useful if X includes 0))
sunflower.ConstPower <- update(sunflower.ols, . ~., weights = varConstPower(form = ~ rescale_ND))
# compare all models by AIC
AIC(sunflower.ols, sunflower.fixed, sunflower.power, sunflower.ConstPower, sunflower.exp)
##                      df      AIC
## sunflower.ols         3 27671.58
## sunflower.fixed       3 34444.82
## sunflower.power       4 27168.59
## sunflower.ConstPower  5 27082.30
## sunflower.exp         4 27077.54
summary(sunflower.exp)
## Generalized least squares fit by REML
##   Model: logSunflower ~ rescale_ND 
##   Data: data 
##        AIC      BIC    logLik
##   27077.54 27108.22 -13534.77
## 
## Variance function:
##  Structure: Exponential of variance covariate
##  Formula: ~rescale_ND 
##  Parameter estimates:
##      expon 
## 0.01781603 
## 
## Coefficients:
##                  Value   Std.Error  t-value p-value
## (Intercept)  0.7471443 0.008183348 91.30057       0
## rescale_ND  -0.0024335 0.000568697 -4.27903       0
## 
##  Correlation: 
##            (Intr)
## rescale_ND -0.841
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -3.1123178 -0.8304163 -0.1021378  0.7680689  3.4945792 
## 
## Residual standard error: 0.4409573 
## Degrees of freedom: 15830 total; 15828 residual